function idat = fsb_normalize_baseline_zscore(idat,sandbox)

% FSB - DEV: Transform data into z-scores
%
% EXAMPLE:
% nidat = fsb_normalize_baseline_zscore(idat,sandbox)
%
% INPUT:
% idat: 4-D image data
% sandbox:      sandbox experiment struct
%
% OUTPUT:
% Diagnostic maps
%
% CALLED BY:
% FSB.m
%
% NOTES:
% work in progress, not fully tested
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
%
%$ Revision 1.0
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

h = waitbar(0,'Dividing trials by own baseline intensity...');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Transform int16 to single for further processing
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

idat = single(idat);
idat(idat==0)=NaN;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Define baseline period and do z-score transform
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
blength=zeros(1,max(sandbox.intrial(:,3)));

vxcorr = zeros(1,max(sandbox.intrial(:,3)));
for x = 1:max(sandbox.intrial(:,3));

    waitbar(x/max(sandbox.intrial(:,3)));
    scanind = find (sandbox.intrial(:,3)==x);
    if min(scanind)==0;
        scanind(1)=[];
    end
    idat_red = idat(:,:,:,scanind); % Gives volumes in trials
    bsl_trial = sum(sandbox.hemodynamics(scanind,1:4),2); % Doing alternative baseline calculation, stimulation regressors only

    % Looks for baseline volumes by looking for values that are below
    % 1/20th of maximum response
    bsl_max = max(bsl_trial);
    bsl_ind = find(bsl_trial<(bsl_max/20));
    %bsl_ind = 1:3; % Use this for constant baseline values

    idat_red = idat_red - repmat(nanmean(idat_red(:,:,:,bsl_ind),4),[1 1 1 size(idat_red,4)]);

    bsl_ind2 = nanstd(idat_red(:,:,:,bsl_ind),1,4);
    blength(x)= size(bsl_ind,1);

    voxel_corr = 1/bsl_ind2; % correction factor voxel by voxel
    voxel_corr_4D = repmat(voxel_corr,[1 1 1 size(idat_red,4)]);
    idat_red = idat_red.* voxel_corr_4D;
    idat(:,:,:,scanind) = idat_red;
    voxel_corr(voxel_corr==inf)=1;
    vxcorr(x) = nanmean(voxel_corr(:));

    disp(['bsl_ind2 : ' num2str(nanmean(bsl_ind2(:)))]);
    disp(['mean voxel_corr : ' num2str(nanmean(voxel_corr(:)))]);
    disp(['min voxel_corr : ' num2str(min(voxel_corr(:)))]);
    disp(['max voxel_corr : ' num2str(max(voxel_corr(:)))]);

end;

disp(['Mean baseline length: ' num2str(mean(blength))]);
disp(['mean voxel corr: ' num2str(mean(vxcorr))]);
close(h);

idat = idat*100;

vxcorr = vxcorr/nanmean(idat(:))*100;
figure; bar(vxcorr);
set(gcf,'Name','Trial zscore intensity correction');
xlabel('Trial');
ylabel('Average signal deviation %');
end
